import pylab as pl #导入绘图模块
from matplotlib.pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']   #显示中文
mpl.rcParams['axes.unicode_minus']=False       #显示负号
import pandas as pd
import numpy as np
from PyEMD import EMD, Visualisation
from dtaidistance import dtw

from scipy.signal import argrelextrema

#进行样条差值
import scipy.interpolate as spi
from scipy import stats

import math
import numpy as np
import matplotlib.pyplot as plt
def moving_average(interval, windowsize):
    window = np.ones(int(windowsize)) / float(windowsize)
    re = np.convolve(interval, window, 'same')
    # re = np.convolve(interval, window, 'full')
    # re = np.convolve(interval, window, 'valid')
    return re
data = pd.read_excel('H033_L1_QDS_DCE_20230823153615675.xlsx')

########################## 136999573300
# coilid = 136999573300
coilid = 137049470100
coilstr = str(coilid)[-4:]

data1 = data[data['coilid'] == coilid]
data1['temp1_2'] = data1['temp1_2'].astype(float)
data11 = data1[data1['temp1_2'] >= 400]
data11 = data11.reset_index(drop=True)
df_max=data11['temp1_2'].max()
df_min=data11['temp1_2'].min()
df0 = data11.copy()
df0['time'] = df0['pos'].astype(float)
df0['temp1_2'] = df0['temp1_2'].astype(float)
timestamp_array = df0['time'].values
value_array = df0['temp1_2'].values
t = timestamp_array
y = value_array
#5个点作滑动平均，平滑曲线
N = 5
#前2本身后2取平均
y_av = moving_average(y, N)
#解决边界点
if N % 2 == 0:
    for i in range(1, int(N / 2) + 1):
        print(i)
        y_av[i-1] = y_av[i-1] * N / (int(N / 2) + i - 1)
    for j in range(1, int(N / 2)):
        print(j)
        y_av[-j] = y_av[-j] * N / (int(N / 2) + j)
    print("偶数")
else:
    for i in range(1, int(N / 2) + 1):
        print(i)
        y_av[i-1] = y_av[i-1] * N / (int(N / 2) + i)
    for j in range(1, int(N / 2) + 1):
        print(j)
        y_av[-j] = y_av[-j] * N / (int(N / 2) + j)
    print("奇数")
df0['new_value'] = y_av



x = [float(v) for v in df0['new_value']]
lx = len(x)
print("共有%d点数据" % (len(x)))
sampling_rate = len(x)  # 取样频率(来自传感器说明书)
fft_size = len(x)  # FFT处理的数据样本数
x = np.asarray(x)
####特征曲线段
# data111 = df0[515:665]
data111 = df0[285:385]


data111['value'] = data111['new_value'].astype(float)
x1 = data111.value.tolist()
l = len(x1)
df_out = pd.DataFrame(columns=['i','pos', 'dtw'])
dict = {}
print(l)
#移动步长
# move_step = 1
move_step = int(l*0.75)
i_max = int((lx-l)/move_step)+1
for i in range(0,i_max):
    print(i)
    data222 = df0[int(i*move_step):int(i*move_step + l)]
    # data222 = df0[i:i+l]
    data222['value'] = data222['new_value'].astype(float)
    x2 = data222.value.tolist()
    distance = dtw.distance(x1, x2)
    dict['i'] = i
    dict['pos'] = int(i*move_step)
    dict['dtw'] = distance
    new_row = pd.Series(dict)
    df_out = df_out.append(new_row, ignore_index=True)
    x11 = data222['pos'].values
    y11 = data222['temp1_2'].values

    x22 = data222['pos'].values
    y22 = data222['value'].values



q1 = df_out['dtw'].quantile(0.25)
q3 = df_out['dtw'].quantile(0.75)
iqr_val = q3 - q1
min3 = q1 - 3 * iqr_val
min15 = q1 - 1.5 * iqr_val
df_out['q1'] = q1
df_out['min3'] = min3
df_out['min15'] = min15

writer = pd.ExcelWriter('dtw_0904_'+coilstr+'_075.xlsx')
df_out.to_excel(writer, sheet_name='Sheet1', index=False)
writer.save()
a1 = stats.percentileofscore(df_out['dtw'],q1)
a2 = stats.percentileofscore(df_out['dtw'],min15)
a3 = stats.percentileofscore(df_out['dtw'],min3)

print(a1)
print(a2)
print(a3)


print(100-a1)
print(a1-a2)
print(a2-a3)
print(a3)


# df_out = df_out[df_out['dtw'] <= min15]
df_out = df_out[df_out['dtw'] <= q1]

df_out = df_out.reset_index(drop=True)
df0['color'] = 'b'
x11 = df0['pos'].values
y11 = df0['temp1_2'].values

plt.figure(figsize=(14, 8))
plt.plot(x11, y11, 'b')
# plt.ylim(550, 650)
y_max=int(df_max*1.1)
y_min=int(df_min*0.9)
plt.ylim(y_min, y_max)

for index, row in df_out.iterrows():
    # print(index)
    # print(row['i'])
    df0.loc[df0.index[int(row['i'] * move_step):int(row['i'] * move_step + l)], 'color'] = 'r'

    data222 = df0[int(row['i']*move_step):int(row['i']*move_step + l)]

    x122 = data222['pos'].values
    y122 = data222['temp1_2'].values
    plt.plot(x122, y122, 'r')

x1122 = data111['pos'].values
y1122 = data111['temp1_2'].values
plt.plot(x1122, y1122, 'black')

plt.xlabel('Pos')
plt.ylabel('Value')
png_name = 'test_0904_'+coilstr+'_075.png'
plt.savefig(png_name)

print('finish')


